//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2012 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : boundaryzetaphi_kernels.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef BOUNDARYZETAPHI_KERNELS_H
#define BOUNDARYZETAPHI_KERNELS_H

namespace wavescomputer
{

namespace boundaryzetaphi
{

namespace reflectiveneumann
{

// Reflective Neumann boundary (North version)
template <class T>
__global__ void kernel_bzp_rn_north(T* phi,
                                    T* zeta,
                                    const T* phiN,
                                    const T* zetaN,
                                    const T* phi_in_N_ext, 
                                    const T* zeta_in_N_ext, 
                                    const T* phi_in_N_cen, 
                                    const T* zeta_in_N_cen, 
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2); // not coalesced

    if (tid < Nx1) 
    {
         phi[tid2D] = phi[tid2D - 1] 
                    + phi_in_N_ext[tid1D] 
                    - phi_in_N_cen[tid1D];

        zeta[tid2D] = zeta[tid2D - 1] 
                    + zeta_in_N_ext[tid1D] 
                    - zeta_in_N_cen[tid1D];
    }
}

// Reflective Neumann boundary (East version)
template <class T>
__global__ void kernel_bzp_rn_east(T* phi,
                                   T* zeta,
                                   const T* phiE,
                                   const T* zetaE,
                                   const T* phi_in_E_ext, 
                                   const T* zeta_in_E_ext, 
                                   const T* phi_in_E_cen, 
                                   const T* zeta_in_E_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2) 
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + Nx1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
         phi[tid2D] = phi[tid2D - ld] 
                    + phi_in_E_ext[tid1D] 
                    - phi_in_E_cen[tid1D];

        zeta[tid2D] = zeta[tid2D - ld] 
                    + zeta_in_E_ext[tid1D] 
                    - zeta_in_E_cen[tid1D];
    }
}

// Reflective Neumann boundary (South version)
template <class T>
__global__ void kernel_bzp_rn_south(T* phi,
                                    T* zeta,
                                    const T* phiS,
                                    const T* zetaS,
                                    const T* phi_in_S_ext, 
                                    const T* zeta_in_S_ext, 
                                    const T* phi_in_S_cen, 
                                    const T* zeta_in_S_cen, 
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH - 1;     // not coalesced

    if (tid < Nx1) 
    {
         phi[tid2D] = phi[tid2D + 1] 
                    + phi_in_S_ext[tid1D] 
                    - phi_in_S_cen[tid1D];

        zeta[tid2D] = zeta[tid2D + 1] 
                    + zeta_in_S_ext[tid1D] 
                    - zeta_in_S_cen[tid1D];
    }
}

// Reflective Neumann boundary (West version)
template <class T>
__global__ void kernel_bzp_rn_west(T* phi, 
                                   T* zeta,
                                   const T* phiW,
                                   const T* zetaW,
                                   const T* phi_in_W_ext, 
                                   const T* zeta_in_W_ext, 
                                   const T* phi_in_W_cen, 
                                   const T* zeta_in_W_cen, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH - 1) * ld + BORDER_WIDTH + tid;     // coalesced

    if (tid < Nx2) 
    {
         phi[tid2D] = phi[tid2D + ld]
                    + phi_in_W_ext[tid1D] 
                    - phi_in_W_cen[tid1D];

        zeta[tid2D] = zeta[tid2D + ld] 
                    + zeta_in_W_ext[tid1D] 
                    - zeta_in_W_cen[tid1D];
    }
}

} // namespace reflectiveneumann




namespace weaklyreflective
{

// weakly-reflective boundary (North version)
template <class T>
__global__ void kernel_bzp_wr_north(T* phi,
                                    T* zeta,
                                    const T* phiN,
                                    const T* zetaN,
                                    const T* phi_in_N_ext,
                                    const T* zeta_in_N_ext,
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2) 
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + (BORDER_WIDTH + Nx2); // not coalesced

    if (tid < Nx1) 
    {
         phi[tid2D] = phiN[tid1D] 
                    + phi_in_N_ext[tid1D];

        zeta[tid2D] = zetaN[tid1D] 
                    + zeta_in_N_ext[tid1D];
    }
}

// weakly-reflective boundary (East version)
template <class T>
__global__ void kernel_bzp_wr_east(T* phi,
                                   T* zeta,
                                   const T* phiE, 
                                   const T* zetaE, 
                                   const T* phi_in_E_ext, 
                                   const T* zeta_in_E_ext, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + Nx1) * ld + BORDER_WIDTH + tid;   // coalesced

    if (tid < Nx2) 
    {
         phi[tid2D] = phiE[tid1D] 
                    + phi_in_E_ext[tid1D];

        zeta[tid2D] = zetaE[tid1D] 
                    + zeta_in_E_ext[tid1D];
    }
}

// weakly-reflective boundary (South version)
template <class T>
__global__ void kernel_bzp_wr_south(T* phi, 
                                    T* zeta,
                                    const T* phiS, 
                                    const T* zetaS, 
                                    const T* phi_in_S_ext, 
                                    const T* zeta_in_S_ext, 
                                    const int ld,
                                    const int Nx1,
                                    const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH + tid) * ld + BORDER_WIDTH - 1;     // not coalesced

    if (tid < Nx1) 
    {
         phi[tid2D] = phiS[tid1D] 
                    + phi_in_S_ext[tid1D];

        zeta[tid2D] = zetaS[tid1D] 
                    + zeta_in_S_ext[tid1D];
    }
}

// weakly-reflective boundary (West version)
template <class T>
__global__ void kernel_bzp_wr_west(T* phi, 
                                   T* zeta,
                                   const T* phiW, 
                                   const T* zetaW, 
                                   const T* phi_in_W_ext, 
                                   const T* zeta_in_W_ext, 
                                   const int ld,
                                   const int Nx1,
                                   const int Nx2)
{
    int tid = bx * Bx + tx;
    int tid1D = tid;                                              // coalesced
    int tid2D = (BORDER_WIDTH - 1) * ld + BORDER_WIDTH + tid;     // coalesced

    if (tid < Nx2) 
    {
         phi[tid2D] = phiW[tid1D] 
                    + phi_in_W_ext[tid1D];

        zeta[tid2D] = zetaW[tid1D] 
                    + zeta_in_W_ext[tid1D];
    }
}

} // namespace weaklyreflective

} // namespace boundaryzetaphi

} // namespace wavescomputer

#endif // BOUNDARYZETAPHI_KERNELS_H

